function [nProd,age,cohort,firmID,quality,ageProd,cohortProd,ProdID,numBS,numCN,check] = simulation_function(p,eq,seed,algS)

rng(seed,'combRecursive');

%% Define Time Period
nSubPer = 1/algS.delt;
nBurn   = algS.nBurn/algS.delt;
nTKeep  = algS.nReal/algS.delt;                      % number of simulation period kept for statistics
nTSim   = nBurn + nTKeep;                            % number of simulation period (total)

%% Preallocation for reporting
nProd    = zeros(algS.maxFirm,algS.nReal*4 + 1);   % firm-specific
age      = zeros(algS.maxFirm,algS.nReal*4 + 1);   % firm-specific
cohort   = zeros(algS.maxFirm,algS.nReal*4 + 1);   % firm-specific
firmID   = zeros(algS.maxFirm,algS.nReal*4 + 1);   % firm-specific
quality  = zeros(algS.maxFirm,algS.NNExt,algS.nReal*4 + 1); % firm-product-specific
ageProd  = zeros(algS.maxFirm,algS.NNExt,algS.nReal*4 + 1); % firm-product-specific
cohortProd  = zeros(algS.maxFirm,algS.NNExt,algS.nReal*4 + 1); % firm-product-specific
ProdID  = zeros(algS.maxFirm,algS.NNExt,algS.nReal*4 + 1); % firm-product-specific
numBS    = zeros(algS.maxFirm,algS.nReal*4 + 1);   % number of business stealing in a year
numCN    = zeros(algS.maxFirm,algS.nReal*4 + 1);   % number of cannibalization in a year


%% Initial Quality Dist

qualityD  = zeros(algS.nFirm,algS.NNExt);
qualityD(:,1) = 1;                             % All firms start with relative quality equal to 1 
ageProdD = -1*ones(algS.nFirm,algS.NNExt);
ageProdD(:,1) = 0;
cohortProdD = -1*ones(algS.nFirm,algS.NNExt);
cohortProdD(:,1) = 0;
ProdIDD = -1*ones(algS.nFirm,algS.NNExt);
ProdIDD(:,1) = (1:algS.nFirm)';

qual_Step = [qualityD(:) ageProdD(:)];
qual_Step = qual_Step((qual_Step(:,2)>-1),:);

nProdD    = ones(algS.nFirm,1);
ageD      = zeros(algS.nFirm,1);
cohortD   = zeros(algS.nFirm,1);
firmIDD   = (1:algS.nFirm)';

maxfirmID = algS.nFirm;
maxprodID = algS.nFirm;
numBSD   = zeros(algS.nFirm,1);   
numCND   = zeros(algS.nFirm,1);   

subPer    = 1; 
period    = 1;

for j = 1:nTSim
    j*algS.delt
    tot_nProdD=sum(nProdD); % total number of product at the beginning
    firmIndex = (1:size(nProdD,1))';
    %% State determination (five mutually exclusive events)
    draw            = rand(size(nProdD,1),1);
    stateBS         = (eq.xE*algS.delt).*nProdD;
    stateCN         = stateBS + (eq.xI*algS.delt).*nProdD;
    stateST         = stateCN + ((eq.xE+eq.xe)*algS.delt).*nProdD;
    stateEX         = stateST + ((eq.xE+eq.xe+eq.xI)*algS.delt).*nProdD;
    if max(stateEX)>=1
        warning('Probability is greater than 1')
    end
    BS = (draw < stateBS);
    CN = logical((1 - BS).*(draw < stateCN));
    ST = logical((1 - BS - CN).*(draw < stateST));
    EX = logical((1 - BS - CN - ST).*(draw < stateEX));
    SQ = logical((1 - BS - CN - ST - EX).*(draw <= 1)); % status quo
    
    nProdD          = min(nProdD + BS + CN - EX,algS.NNExt - 1);      % update number of product line    
    ageD            = ageD + 1;
    
    %% Update qualities
    
    % business stealing
    newQualBS       = datasample(qual_Step,sum(BS));
    IndBS           = (nProdD(BS) - 1)*size(nProdD,1) + firmIndex(BS);
    qualityD(IndBS) = newQualBS(:,1) * (1+p.lambdaE);
    cohortProdD(IndBS) = j*ones(size(newQualBS,1),1);
    ProdIDD(IndBS) = (1:size(newQualBS,1))'+maxprodID;
    maxprodID=maxprodID+size(newQualBS,1);
    
    firmIndBS = firmIndex(BS);
    ageProdD(BS,:)=(ageProdD(BS,:)>-1)+ageProdD(BS,:); % prod age updates
    ageProdD(IndBS) = 0;
    
    % cannibalization
    drawCN  = rand(sum(CN),1);       % we randomly drop the quality
    dropQualsCN = ceil(drawCN.*(nProdD(CN)-1));
    firmIndCN = firmIndex(CN);
    ageProdD(CN,:)=(ageProdD(CN,:)>-1)+ageProdD(CN,:); % prod age updates
    % step 1: new better product
    oldQual=zeros(sum(CN),1);
    for s=1:sum(CN)
        oldQual(s,1)=qualityD(firmIndCN(s),dropQualsCN(s));
    end
    IndCN           = (nProdD(CN) - 1)*size(nProdD,1) + firmIndex(CN);
    qualityD(IndCN) = oldQual(:,1) * (1+p.lambdaI);
    cohortProdD(IndCN) = j*ones(size(oldQual,1),1);
    ProdIDD(IndCN) = (1:size(oldQual,1))'+maxprodID;
    maxprodID=maxprodID+size(oldQual,1);
    ageProdD(IndCN) = 0;
    % step 2: old discounted product
    for s = 1:sum(CN)
        qualityD(firmIndCN(s),dropQualsCN(s))=(1-p.lambdaC)*qualityD(firmIndCN(s),dropQualsCN(s));
    end
    
    % being stolen
    drawST  = rand(sum(ST),1);       % we randomly drop the quality
    dropQualsST = ceil(drawST.*(nProdD(ST)));
    firmIndST = firmIndex(ST);
    for s = 1:sum(ST)
        qualityD(firmIndST(s),dropQualsST(s))=(1-p.lambdaS)*qualityD(firmIndST(s),dropQualsST(s));
    end
    ageProdD(ST,:)=(ageProdD(ST,:)>-1)+ageProdD(ST,:); % prod age updates
    
    % status quo
    firmIndSQ = firmIndex(SQ);
    ageProdD(SQ,:)=(ageProdD(SQ,:)>-1)+ageProdD(SQ,:); % prod age updates
    
    % product exit
    drawEX  = rand(sum(EX),1);       % we randomly drop the quality
    dropQualsEX = ceil(drawEX.*(nProdD(EX)+1));
    firmIndEX = firmIndex(EX);
    ageProdD(EX,:)=(ageProdD(EX,:)>-1)+ageProdD(EX,:); % prod age updates
    for s = 1:sum(EX)
        qualityD(firmIndEX(s),dropQualsEX(s):end-1) = qualityD(firmIndEX(s),dropQualsEX(s)+1:end);  % shifting such that one randomly quality drops (this is kind of problematic if all quality matrix fills with product lines)
        ageProdD(firmIndEX(s),dropQualsEX(s):end-1) = ageProdD(firmIndEX(s),dropQualsEX(s)+1:end);
        cohortProdD(firmIndEX(s),dropQualsEX(s):end-1) = cohortProdD(firmIndEX(s),dropQualsEX(s)+1:end);
        ProdIDD(firmIndEX(s),dropQualsEX(s):end-1) = ProdIDD(firmIndEX(s),dropQualsEX(s)+1:end);
    end
    
    % count number of BS and CN
    numBSD(BS)=numBSD(BS)+1;
    numCND(CN)=numCND(CN)+1;
    
    % Check if firm exits
    ifExit           = logical((nProdD<1).*EX);
    firmIndExit      = firmIndex(ifExit);
    for s=1:sum(ifExit)
        nProdD(firmIndExit(s))=NaN;
        ageD(firmIndExit(s))=NaN;
        cohortD(firmIndExit(s))=NaN;
        firmIDD(firmIndExit(s))=NaN;
        qualityD(firmIndExit(s),:)=NaN;
        ageProdD(firmIndExit(s),:)=NaN;
        cohortProdD(firmIndExit(s),:)=NaN;
        ProdIDD(firmIndExit(s),:)=NaN;
        numBSD(firmIndExit(s),:)=NaN;
        numCND(firmIndExit(s),:)=NaN;
    end
    nProdD=rmmissing(nProdD); % remove missing entries
    ageD=rmmissing(ageD); % remove missing entries
    cohortD=rmmissing(cohortD); % remove missing entries
    firmIDD=rmmissing(firmIDD); % remove missing entries
    qualityD=rmmissing(qualityD); % remove missing entries
    ageProdD=rmmissing(ageProdD); % remove missing entries
    cohortProdD=rmmissing(cohortProdD); % remove missing entries
    ProdIDD=rmmissing(ProdIDD); % remove missing entries
    numBSD=rmmissing(numBSD); % remove missing entries
    numCND=rmmissing(numCND); % remove missing entries
    
    % product/firm entry
    % i. update firm variables
    NumNewP=ceil(eq.xe*algS.delt*tot_nProdD);
    nProdD=[nProdD;ones(NumNewP,1)];
    ageD=[ageD;zeros(NumNewP,1)];
    cohortD=[cohortD;j*ones(NumNewP,1)];
    firmIDD=[firmIDD;(1:NumNewP)'+maxfirmID];
    maxfirmID=maxfirmID+NumNewP; % track the lastest firm ID
    numBSD=[numBSD;zeros(NumNewP,1)];
    numCND=[numCND;zeros(NumNewP,1)];
    % ii. update product variables
    NewqualityD  = zeros(NumNewP,algS.NNExt);
    newQualEN       = datasample(qual_Step,NumNewP);
    %newQualEN       = ones(NumNewP,1);
    NewqualityD(:,1)  = newQualEN(:,1) * (1+p.lambdae);
    NewageProdD  = -1*ones(NumNewP,algS.NNExt);
    NewageProdD(:,1)  = zeros(NumNewP,1);
    qualityD=[qualityD;NewqualityD];
    ageProdD=[ageProdD;NewageProdD];
    NewcohortProdD = -1*ones(NumNewP,algS.NNExt);
    NewcohortProdD(:,1)  = j*ones(NumNewP,1);
    cohortProdD = [cohortProdD;NewcohortProdD];
    NewProdIDD = -1*ones(NumNewP,algS.NNExt);
    NewProdIDD(:,1)  = (1:NumNewP)'+maxprodID;
    maxprodID=maxprodID+NumNewP;
    ProdIDD = [ProdIDD;NewProdIDD];
    
    % New quality distribution
    qual_Step = [qualityD(:) ageProdD(:)];
    qual_Step = qual_Step((qual_Step(:,2)>-1),:);
    check(j,1)=mean(qual_Step(:,1)); % report q_bar
    check(j,2)=sum(qual_Step(:,1).^5)^(1/5); % report Q
    check(j,3)=max(nProdD); % maximum number of products per firm
    check(j,4)=size(nProdD,1); % number of firm
    check(j,5)=sum(nProdD==algS.NNExt-1)/size(nProdD,1)*100; % share of firms at the upper bound (%)

    %% Save Outcome
    if j==nBurn % this is intitial distribution
        nProd(1:size(nProdD,1),period)     = nProdD; % firm-specific
        age(1:size(nProdD,1),period)       = ageD;   % firm-specific
        cohort(1:size(nProdD,1),period)       = cohortD;   % firm-specific
        firmID(1:size(nProdD,1),period)       = firmIDD;   % firm-specific
        quality(1:size(nProdD,1),:,period) = qualityD; % firm-product-specific
        ageProd(1:size(nProdD,1),:,period) = ageProdD; % firm-product-specific
        cohortProd(1:size(nProdD,1),:,period) = cohortProdD; % firm-product-specific
        ProdID(1:size(nProdD,1),:,period) = ProdIDD; % firm-product-specific
        period              = period + 1;
        numBSD   = zeros(size(nProdD,1),1);   
        numCND   = zeros(size(nProdD,1),1); 
    end
    
    if j>nBurn % keeping period

    if subPer == 1/4*nSubPer % end of quarter 1
            % Reporting the end-quarter values
            nProd(1:size(nProdD,1),period)         = nProdD; 
            age(1:size(nProdD,1),period)           = ageD;
            cohort(1:size(nProdD,1),period)           = cohortD;
            firmID(1:size(nProdD,1),period)           = firmIDD;
            quality(1:size(nProdD,1),:,period)     = qualityD;
            ageProd(1:size(nProdD,1),:,period)     = ageProdD;
            cohortProd(1:size(nProdD,1),:,period)     = cohortProdD;
            ProdID(1:size(nProdD,1),:,period)     = ProdIDD;
            numBS(1:size(nProdD,1),period)   = numBSD;
            numCN(1:size(nProdD,1),period)   = numCND;
            % reset
            numBSD   = zeros(size(nProdD,1),1);   
            numCND   = zeros(size(nProdD,1),1);   
    end
    if subPer == 2/4*nSubPer % end of quarter 2
            % Reporting the end-quarter values
            nProd(1:size(nProdD,1),period+1)         = nProdD; 
            age(1:size(nProdD,1),period+1)           = ageD;
            cohort(1:size(nProdD,1),period+1)           = cohortD;
            firmID(1:size(nProdD,1),period+1)           = firmIDD;
            quality(1:size(nProdD,1),:,period+1)     = qualityD;
            ageProd(1:size(nProdD,1),:,period+1)     = ageProdD;
            cohortProd(1:size(nProdD,1),:,period+1)     = cohortProdD;
            ProdID(1:size(nProdD,1),:,period+1)     = ProdIDD;
            numBS(1:size(nProdD,1),period+1)   = numBSD;
            numCN(1:size(nProdD,1),period+1)   = numCND;
            % reset
            numBSD   = zeros(size(nProdD,1),1);   
            numCND   = zeros(size(nProdD,1),1);   
    end
    if subPer == 3/4*nSubPer % end of quarter 3
            % Reporting the end-quarter values
            nProd(1:size(nProdD,1),period+2)         = nProdD; 
            age(1:size(nProdD,1),period+2)           = ageD;
            cohort(1:size(nProdD,1),period+2)           = cohortD;
            firmID(1:size(nProdD,1),period+2)           = firmIDD;
            quality(1:size(nProdD,1),:,period+2)     = qualityD;
            ageProd(1:size(nProdD,1),:,period+2)     = ageProdD;
            cohortProd(1:size(nProdD,1),:,period+2)     = cohortProdD;
            ProdID(1:size(nProdD,1),:,period+2)     = ProdIDD;
            numBS(1:size(nProdD,1),period+2)   = numBSD;
            numCN(1:size(nProdD,1),period+2)   = numCND;
            % reset
            numBSD   = zeros(size(nProdD,1),1);   
            numCND   = zeros(size(nProdD,1),1);   
    end
    if subPer == nSubPer % end of quarter 4
            % Reporting the end-quarter values
            nProd(1:size(nProdD,1),period+3)         = nProdD; 
            age(1:size(nProdD,1),period+3)           = ageD;
            cohort(1:size(nProdD,1),period+3)           = cohortD;
            firmID(1:size(nProdD,1),period+3)           = firmIDD;
            quality(1:size(nProdD,1),:,period+3)     = qualityD;
            ageProd(1:size(nProdD,1),:,period+3)     = ageProdD;
            cohortProd(1:size(nProdD,1),:,period+3)     = cohortProdD;
            ProdID(1:size(nProdD,1),:,period+3)     = ProdIDD;
            numBS(1:size(nProdD,1),period+3)   = numBSD;
            numCN(1:size(nProdD,1),period+3)   = numCND;
            % reset
            numBSD   = zeros(size(nProdD,1),1);   
            numCND   = zeros(size(nProdD,1),1);   
            period = period + 4;
            subPer = 1;
        else
            subPer = subPer + 1;
    end
end
end
end

